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Abstract 

We present an investigation of the dependence of searches for boosted Higgs bosons 
using jet substructure on the perturbative and non-perturbative parameters of the 
Herwig++ Monte Carlo event generator. Values are presented for a new tune of the 
parameters of the event generator, together with the an estimate of the uncertainties 
based on varying the parameters around the best-fit values. 



1 Introduction 

Monte Carlo simulations are an essential tool in the analysis of modern collider 
experiments. These event generators contain a large number of both perturbative 
and non-perturbative parameters which are tuned to a wide range of experimental 
data. While significant effort has been devoted to the tuning of the parameters to 
produce a best fit there has been much less effort understanding the uncertainties 
in these results. Historically a best fit result, or at best a small number of tunes, 
are produced and used to predict observables making it difficult to assess the 
uncertainty on any prediction. The "Perugia" tunes [TJ|2] have addressed this 
by producing a range of tunes by varying specific parameters in the Pythia [3] 
event generator to produce an uncertainty. 

Here we make use of the Professor Monte Carlo tuning system [2] to give 
an assessment of the uncertainty by varying all the parameters simultaneously 
about the best-fit values by diagonalizing the error matrix. This then allows us 
to systematically estimate the uncertainty on any Monte Carlo prediction from 
the tuning of the event generator. We will illustrate this by considering the 
uncertainty on jet substructure searches for the Higgs boson at the LHC. 

As the LHC takes increasing amounts of data the discovery of the Higgs boson 
is likely in the near future. Once we have discovered the Higgs boson, most likely 
in the diphoton channel, it will be vital to explore other channels and determine 
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if the properties of the observed Higgs boson are consistent with the Standard 
Model. For many years it was believed that it would be difficult, if not impossible, 
to observe the dominant h° — >■ bb decay mode of a light Higgs boson. However, 
in recent years the use of jet substructure (5 -20 offers the possibility of observing 
this mode. Jet substructure for h° — > bb as a Higgs boson search channel, was 
first studied in Ref. |5] building on previous work of a heavy Higgs boson decaying 



to W ± bosons 16 
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and SUSY decay chains 
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high-energy WW scattering 
and subsequently reexamined in Refs. (8 , 15 . Recent studies at the LHC 23 
have also shown this approach to be promising. 

The study in Ref. was carried out using the (FORTRAN) HERWIG 6.510 
event generator 26 27 together with the simulation of the underlying event using 
JIMMY 4.31 
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In order to allow the inclusion of new theoretical developments 
and improvements in non-perturbative modelling a new simulation based on the 
same physics philosophy Herwig++, currently version 2.6 |29|[30] , is now preferred 
for the simulation of hadron-hadron collisions. 

Herwig++ includes both an improved theoretical description of perturbative 
QCD radiation, in particular for radiation from heavy quarks, such as bottom, 
together with improved non-perturbative modeling, especially of multiple parton- 
parton scattering and the underlying event. In FORTRAN HERWIG a crude 
implementation of the dead-cone effect 31 meant that there was no radiation 



from heavy quarks for evolution scales below the quark mass, rather than a 
smooth suppression of soft collinear radiation. In Herwig++ an improved choice 
of evolution variable [32] allows evolution down to zero transverse momentum for 
radiation from heavy particles and reproduces the correct soft limit. There have 
also been significant developments of the multiple-parton scattering model of the 
underlying event 33 34 , including colour reconnections 35 and tuning to LHC 
data 
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The background to jet substructure searches for the Higgs boson comes from 
QCD jets which mimic the decay of a boosted heavy particle. Although Her- 
wig++ has performed well in some early studies of jet substructure |25|[37][38] , 
it is important that we understand the uncertainties in our modelling of the 
background jets which lie at the tail of the jet mass distribution. 

In addition we improve the simulation of Higgs boson decay by implement- 
ing the next-to-leading-order (NLO) corrections to Higgs boson decay to heavy 
quarks in the POWHEG [39||40] formalism. 



In the next section we present our approach for the tuning of the parameters, 
which effect QCD radiation and hadronization, in Herwig+-|- together with the 
results of our new tune. We then recap the key features of the Butterworth, 
Davison, Rubin and Salam (BDRS) jet substructure technique of Ref. (5]. This 
is followed by our results using both the leading and next-to-leading-order matrix 
elements in Herwig++ with implementation of the next-to-leading-order Higgs 
boson decays and our estimate on the uncertainties. 



2 Tuning Herwig++ 



3 



2 Tuning Herwig++ 

Any jet substructure analysis is sensitive to changes in the simulation of initial- 
and final-state radiation, and hadronization. In particular the non-perturbative 
nature of the phenomenological hadronization model means there are a number of 
parameters which are tuned to experimental results. Herwig++ uses an improved 



angular-ordered parton shower algorithm 29 32 to describe perturbative QCD 
radiation together with a cluster hadronization model 29 41 . 

The Herwig++ cluster model is based on the concept of preconfinement 42 . 
At the end of the parton-shower evolution all gluons are non-perturbatively split 
into quark-antiquark pairs. All the partons can then be formed into colour- 
singlet clusters which are assumed to be hadron precursors and decay according 
to phase space into the observed hadrons. There is a small fraction of heavy 
clusters for which this is not a reasonable approximation which are therefore first 
fissioned into lighter clusters. The main advantage of this model, when coupled 
with the angular-ordered parton shower is that it has fewer parameters than the 

3] event generator yet still gives a 
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string model as implemented in the Pythia 
reasonable description of collider observables 

To tune Herwig++, and investigate the dependency of observables on the 
shower and hadronization parameters, the Professor Monte Carlo tuning sys- 
tem [4] was used. Professor uses the Rivet analysis framework (44] and a number 
of simulated event samples, with different Monte Carlo parameters, to parame- 
terise the dependence of each observable^ used in the tuning on the parameters 
of the Monte Carlo event generator. A heuristic chi-squared function 

x (p) = 2^wo2^ ^2 > C 1 ) 



o beo 



is constructed where p is the set of parameters being tuned, O are the observables 
used each with weight wq, b are the different bins in each observable distribu- 
tion with associated experimental measurement TZb, error Aj, and Monte Carlo 
prediction f b (p)- Weighting of those observables for which a good description of 
the experimental result is important is used in most cases. The parameterisation 
of the event generator response, f{p), is then used to minimize the x' 2 an d find 
the optimum parameter values. 

There are ten main free parameters which affect the shower and hadronization 
in Herwig+-|-. These are shown in Table [T] along with their default values and 
allowed ranges. 

The gluon mass, GluonMass, is required to allow the non-perturbative decay of 
gluons into qq pairs and controls the energy release in this process. PSplitLight, 



1 Normally this is either an observation such as a multiplicity or a bin in a measured distri- 
bution. 
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Parameter 


Default Value 


Allowed Range 


Scanned Range 


Optimum Value 


Qmin 


0.935 


> o 


0.500- 2.500 


Fixed at default 


GluonMass 


0.95 


0-1 


0.75-1.00 


Fixed at default 


CISmrLight 


0.78 


0-2 


0.30-3.00 


Fixed at default 


CIPowLight 


1.28 


0-10 


0.50-4.00 


Fixed at default 


pTmin 


1.00 


> o 


0.50- 1.50 


0.88 


AlphaMZ 


0.12 


> o 


0.10-0.12 


0.11 


CIMaxLight 


3.25 


0-10 


3.00-4.20 


3.60 


PSplitLight 


1.20 


0-10 


1.00-2.00 


0.90 


PwtDIquark 


0.49 


0-10 


0.10-0.50 


0.33 


PwtSquark 


0.68 


0-10 


0.50-0.80 


0.64 



Tab. 1: The ten parameters to which the jet substructure is most sensitive with 
their default values, the allowed range of these values in Herwig++, the 
range scanned over and the new optimum value found from minimizing 

x' 2 - 



CIPowLight and CIMaxLight control the mass distributions of the clusters pro- 
duced during the fission of heavy clusters. CISmrLight controls the smearing of 
the direction of hadrons containing a (anti) quark from the perturbative evolu- 
tion about the direction of the (anti) quark. AlphaMZ is strong coupling at the 
Z° boson mass and controls the amount of QCD radiation in the parton shower, 
while Qmin controls the infrared behaviour of the strong coupling. pTmin is the 
minimum allowed transverse momentum in the parton shower and controls the 
amount of radiation and the scale at which the perturbative evolution termi- 
nates. PwtDIquark and PwtSquark are the probabilities of selecting a diquark- 
antidiquark or ss quark pair from the vacuum during cluster splitting, and affect 
the production of baryons and strange hadrons respectively. 

Previous experience of tuning Herwig++ has found that Qmin, GluonMass, 
CISmrLight and CIPowLight to be flat, and so it was chosen to fix these at their 



default values 29 



To determine the allowed variation of these parameters Professor was used to 
tune the variables in Table [T] to the observables and weights found in Appendix [A] 
in Tables [5j [6j [7]and[8j The dependence of x' 2 011 the various parameters, about 
the minimum x' 2 value, is then diagonalized. 

The variation of the parameters along the eigenvectors in parameter space 
obtained corresponding to a certain change, Ax' 2 , in x' 2 can then be used to 
predict the uncertainty in the Monte Carlo predictions for specific observables. 

In theory, if the x' 2 measure for the parameterised generator response is ac- 
tually distributed as a true x 2 ■, then a change in the goodness of fit of one will 
correspond to a one sigma deviation from the minima, i.e. the best tune. In 
practice, even the best tune does not fit the data ideally and nor is the x' 2 mea- 
sure actually distributed according to a true x 2 distribution. This means that 
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Fig. 1: The x'V-^df distributions for the parameters that were varied from their 
default values whilst determining the error tune. The scatter of the results 
gives a representation of the systematics of tuning procedure. 

one cannot just use Professor to vary the parameters about the minima to a given 
deviation in the x' 2 measure without using some subjective opinion on the quality 
of the results. 

We simulated one thousand event samples with different randomly selected 
values of the parameters we were tuning. Six hundred of these were used to 
interpolate the generator response. All the event samples were used to select two 
hundred samples randomly two hundred times in order to assess the systematics 
of the interpolation and tuning procedure. A cubic interpolation of the generator 
response was used as this has been shown to give a good description of the Monte 
Carlo behaviour in the region of best generator response R] . The parameters were 
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varied between values shown in Table [TJ The quality of the interpolation was 
checked by comparing the x'V-^dfj where iVdf is the number of observable bins 
used in the tune, in the allowed parameter range on a parameter by parameter 
basis for the observables by comparing the interpolation response with actual 
generator response at the simulated parameter values. Bad regions were removed 
and the interpolation repeated leaving a volume in the 5-dimensional parameter 
space where the interpolation worked well. 

Fig. [I] shows the x' 2 /^df distributions for two hundred tunes based on two 
hundred randomly selected event samples points for the cubic interpolation. The 
spread of these values gives an idea of the systematics of the tuning process 
showing that we have obtained a good fit for our parameterisation of the generator 
response. 

The line indicates the tune which is based on a cubic interpolation from six 
hundred event samples. It is this interpolation which was used to vary x' 2 about 
the minimum to assess the uncertainty on the measured distributions. During 
the tune it was discovered that PSplitLight was relatively insensitive to the 
observables used in the tune. As such, PSplitLight was fixed at the default 
value of 1.20 during the tune and subsequent x' 2 variation. 

Professor was used to vary x' 2 about the minimum value, as described above, 
determining the allowed range for the parameters. As five parameters were even- 
tually varied, there are 10 new sample points - one for each of the parameters 
and one "+" and one "-" along each eigenvector direction in parameter space. 

Various changes in x' 2 were examined. A variation of Ax' 2 /N^ = 5 keeps 
the LEP data within reasonable limits while a variation of A% /2 /A^jf = 10 is too 
large. The values for both Ax' 2 /A^f = 5 and Ax' 2 /N^ = 10 are shown in are 
shown in Tables [2] and Tables [3] respectively. 

The Professor tune was then compared with the internal Herwig++ tuning 
procedure [29] as not all analyses that are in the internal Herwig++ tuning system 
are available in Rivet and subsequently accessible to Professor. Looking at Fig. [2] 
it is found that PSplitLight at a value of 0.90 is favoured and gives a significant 
reduction in the x' 2 /Adf- It was therefore decided to use the values obtained from 
minimisation procedure, but using the value of 0.90 for PSplitLight to maintain 
a good overall description of the data. The new minima for the QCD parameters 
are summarized in the Table [TJ Examples of the new tune and the uncertainty 
band are shown in Figs. [2] and [3] for the out-of-plane transverse momentum and 
thrust measured by DELPHI |45| . 

These error tune values can now be used to predict the uncertainty from 
the tuning of the shower parameters on any observable. In the next section we 
will present an example of using these tunes to estimate the uncertainty on the 
predictions for searches for the Higgs boson using the BDRS jet substructure 
method. 
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Fig. 2: Results from the DELPHI 45 



analysis of out-of-plane pr with-respect- 
to the thrust axis and 1-thrust showing the new tune and the envelopes 
corresponding to a change in Ax' 2 /N^ = 5. 
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Fig. 3: Results from the DELPHI |45j analysis of out of plane pt with- respect- 
to the thrust axis and 1-thrust showing the new tune and the envelopes 
corresponding to a change in Ax' 2 /iV*df = 10. 

3 Jet Substructure Boosted Higgs 

The analysis of Ref. [5] uses a number of different channels for the production of 
the Higgs boson decaying to bb in association with an electroweak gauge boson, 
i.e. the production of h°Z° and h°W ± . Ref. [5] uses the fact that the Higgs boson 
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Parameter 


Direction 


1 


2 


3 


4 


5 


+ 




+ 




+ 




+ 




+ 




pTmin 


0.88 


0.88 


0.88 


0.88 


0.84 


0.93 


0.87 


0.90 


0.89 


0.87 


AlphaMZ 


0.11 


0.11 


0.10 


0.12 


0.12 


0.11 


0.12 


0.11 


0.12 


0.11 


CIMaxLight 


3.61 


3.61 


3.61 


3.61 


3.60 


3.62 


3.66 


3.55 


3.54 


3.67 


PwtDIquark 


0.46 


0.23 


0.33 


0.33 


0.33 


0.33 


0.33 


0.33 


0.33 


0.33 


PwtSquark 


0.64 


0.64 


0.64 


0.64 


0.64 


0.64 


0.62 


0.67 


0.51 


0.78 



Tab. 2: The five directions corresponding to the error tune for a A% /2 /A^df = 5 
and the values the parameters take in each direction. 



Parameter 


Direction 


1 


2 


3 


4 


5 


+ 




+ 




+ 




+ 




+ 




pTmin 


0.88 


0.88 


0.88 


0.88 


0.82 


0.95 


0.86 


0.90 


0.89 


0.87 


AlphaMZ 


0.11 


0.11 


0.10 


0.12 


0.12 


0.10 


0.12 


0.10 


0.12 


0.11 


CIMaxLight 


3.61 


3.61 


3.61 


3.61 


3.59 


3.63 


3.68 


3.52 


3.52 


3.70 


PwtDIquark 


0.51 


0.19 


0.33 


0.33 


0.33 


0.33 


0.33 


0.33 


0.33 


0.33 


PwtSquark 


0.64 


0.64 


0.64 


0.64 


0.65 


0.64 


0.61 


0.68 


0.46 


0.84 



Tab. 3: The five directions corresponding to the error tune for a Ax'V^df = 10 
and the values the parameters take in each direction. 




Fig. 4: A scan of PSplitLight using the internal Herwig++ tuning system with 
the other parameters fixed at their new tuned value. From the total 
X /2 /Ndf we see that a value of 0.90 for PSplitLight is favoured at the 
new tuned parameters driven by the multiplicities. 
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Fig. 5: Results for the reconstructed Higgs boson mass distribution using leading- 
order matrix elements. A SM Higgs boson was assumed with a mass of 
115 GeV. In addition to the full result the contribution from top quark 
pair production (tt), the production of W ± (W+Jet) and Z° (Z+Jet) 
bosons in association with a hard jet, vector boson pair production (VV) 
and the production of a vector boson in association with the Higgs bo- 
son (V+Higgs), are shown. 



predominantly decays to bb in a jet substructure analysis to extract the signal of 
a boosted Higgs boson above the various backgrounds. Their study found that 
the Cambridge-Aachen algorithm 46,47 with radius parameter R = 1.2 gave 
the best results when combined with their jet substructure technique. For our 
study, we used the Cambridge-Aachen algorithm as implemented in the FastJet 
package 48 . Three different event selection criteria are used: 



(a) a lepton pair with 80 GeV < mi+i- < 100 GeV and pr > Pr m to select events 
for Z° 



(b) missing transverse momentum fl > p^ m to select events with Z° — > uu; 
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Fig. 6: Results for the reconstructed Higgs boson mass distribution using leading- 
order matrix elements for top quark pair production (tt), and the produc- 
tion of W ± (W+Jet) and Z° (Z+Jet) bosons in association with a hard 
jet. The next-to-leading-order corrections are included for vector boson 
pair production (VV) and the production of a vector boson in associa- 
tion with the Higgs boson (V+Higgs) as well as in the decay of the Higgs 
boson, h° bb. A SM Higgs boson was assumed with a mass of 115 GeV. 

(c) missing transverse momentum p > 30 GeV and a lepton with pj- > 30 GeV 
consistent with the presence of a W boson with px > p™ in to select events 
with W -4 iv\ 

where pf n = 200 GeV. 

In addition the presence of a hard jet with p T . > p™ m with substructure is 
required. The substructure analysis of Ref. |5] proceeds with the hard jet j with 
some radius Rj, a mass rrij and in a mass-drop algorithm: 

1. the two subjets which were merged to form the jet, ordered such that the 
mass of the first jet rrij 1 is greater than that of the second jet 
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Significance 


Process 


Order 


s 

VB 

Herwig++ default 


s 

VB 

Herwig++ tune 


r~7C\ i —J t 

Z° ->■ l + l 


LO 
NLO 


1.17 
1.57 


1.241°,?? 
l-96l°|° 


Z° ->■ vv 


LO 
NLO 


2.18 
2.95 


2.89l°i° 
4-041°;^ 


W -+lv 


LO 
NLO 


1.88 
2.63 


9 09+O.15 
Z ' OZ -0.27 

3 20 + °- 29 

°- zu -0.36 


Total 


LO 
NLO 


2.98 
4.09 


71 +0.29 
°- ' ± -0.53 

5 20+°- 43 

°- zu -0.81 



Tab. 4: The significance of the different processes for the leading- and next-to- 
leading-order matrix elements. The significance is calculated using all 
masses in the range 112-120 GeV. 



obtained; 
2. if < jirrij and 



mm(p 2 Tji , p 2 Tj2 ) 

—2 AK h,h > J/cut, 



y = ^ in A<, 2 > ycut, (2) 

where AR 2 ji h = {y h - y h ) 2 + - j2 ) 2 , and p Tjl 2 , r] jl 2 , (p jl 2 are the 
transverse momenta, rapidities and azimuthal angles of jets 1 and 2, re- 
spectively, then j is in the heavy particle region. If the jet is not in the 
heavy particle region the procedure is repeated using the first jet. 

This algorithm requires that are 6-tagged and takes \i = 0.67 and y cut = 0.09. 
A uniform 6-tagging efficiency of 60% was used with a uniform mistagging prob- 
ability of 2%. The heavy jet selected by this procedure is considered to be the 
Higgs boson candidate jet. Finally, there is a filtering procedure on the Higgs 
boson candidate jet, j. The jet, j, is resolved on a finer scale by setting a 
new radius i?mt = min(0.3, R b i/2), where from the previous mass-drop condition, 
R b l = AR 2 i j2 . The three hardest subjects of this filtering process are taken to 
be the Higgs boson decay products, where the two hardest are required to be 
6-tagged. 

All three analyses require that: 

• after the reconstruction of the vector boson, there are no additional leptons 
with pseudorapidity \rj\ < 2.5 and px > 30 GeV; 
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Fig. 7: Results for the reconstructed Higgs boson mass distribution using leading- 
order matrix elements. A SM Higgs boson was assumed with a mass of 
115 GeV. The envelope shows the uncertainty from the Monte Carlo 
simulation. 



• other than the Higgs boson candidate, there are no additional 6-tagged jets 
with pseudorapidity \r)\ < 2.5 and px > 50 GeV. 

In addition, due to top contamination, criterion Q requires that other than 
the Higgs boson candidate, there are no additional 6-tagged jets with \t)\ < 3 
and px > 30 GeV. For all events, the candidate Higgs boson jet should have 
The analyses were implemented using the Rivet system 
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Pt > Pt 

The plots shown in Fig. [5] use the leading-order matrix elements for the pro- 
duction and decay of Higgs boson but the W, Z and top [49] have matrix element 
corrections for their decays. The plots shown in Fig. [6] have leading-order ti pro- 
duction, leading-order vector boson plus jet production (with the same matrix 
element corrections as the LO matrix elements) but the NLO vector boson pair 
production 



50 



51 



and NLO vector and Higgs boson associated production 
In addition we have implemented the corrections to the decay h° — > bb in the 
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Fig. 8: Results for the reconstructed Higgs boson mass distribution using leading- 
order matrix elements for top quark pair production, and the production 
of W ± and Z° bosons in association with a hard jet. The next-to-leading- 
order corrections are included for vector boson pair production and the 
production of a vector boson in association with the Higgs boson as well 
as in the decay of the Higgs boson, h° — > bb. A SM Higgs boson was 
assumed with a mass of 115 GeV. The envelope shows the uncertainty 
from the Monte Carlo simulation. 



POWHEG scheme, as described in Appendix |B} The signal significances are out- 
lined in Table HI 

The uncertainties due to the Monte Carlo simulation are shown as bands in 
Figs. [7] and [8| As there are correlations between the different processes the uncer- 
tainty is determined for the sum of all processes. Whilst it would be possible to 
show the envelope for the individual processes, this would not offer any informa- 
tion on the envelope for the sum of the processes which is the result of interest. 
In addition the uncertainty on the significance is shown in Table |4j 
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4 Conclusions 

Monte Carlo simulations are an essential tool in the analysis of modern collider 
experiments. While significant effort has been devoted to the tuning of the pa- 
rameters to produce a best fit there has been much less effort understanding the 
uncertainties in these results. In this paper we have produced a set of tunes which 
can be used to assess this uncertainty using the Herwig++ Monte Carlo event 
generator. 

We then used these tunes to assess the uncertainties on the mass-drop analysis 
of Ref. [5] using Herwig++ with both leading- and next-to-leading-order matrix 
elements including a POWHEG simulation of the decay h° — > bb. 

We find that while the jet substructure technique has significant potential as a 
Higgs boson discovery channel, we need to be confident of our tunes to investigate 
this with Monte Carlo simulations. 

The error tunes and procedure here can now be used in other analyses where 
the uncertainty due to the Monte Carlo simulation is important. 
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A Observables and weights used to tune Herwig++ 

The weights and observables used in the Professor tuning system are outlined in 
the tables below. 



Observable 


Weight 


Observable 


Weight 


ir >± (892) spectrum 


1.0 


A° spectrum 


1.0 


p spectrum 


1.0 


7T spectrum 


1.0 


a; (782) spectrum 


1.0 


p spectrum 


1.0 


S~ spectrum 


1.0 


77' spectrum 


1.0 


if* (892) spectrum 


1.0 


S°(1530) spectrum 


1.0 


(j) spectrum 


1.0 


tt ± spectrum 


1.0 


£±(1385) spectrum 


1.0 


r\ spectrum 


1.0 


7 spectrum 


1.0 


K° spectrum 


1.0 


K spectrum 


1.0 







Tab. 5: Observables used in the tuning and associated weights for observables 
taken from 52 . 
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Observable 


Weight 


Sphericity, S 


1.0 


Energy-energy correlation, EEC 


1.0 


Aplanarity, A 


2.0 


Mean out-of-plane p± in GeV w.r.t. thrust axes vs. x p 


1.0 


Mean charged multiplicity 


150.0 


Mean p± in GeV vs. x p 


1.0 


Planarity, P 


1.0 


Thrust major, M 


1.0 


Oblateness = M — m 


1.0 


Out-of-plane p± in GeV w.r.t. sphericity axes 


1.0 


D parameter 


1.0 


1 - Thrust 


1.0 


Out-of-plane p± in GeV w.r.t. thrust axes 


1.0 


Log of scaled momentum, log(l/x p ) 


1.0 


In-plane p± in GeV w.r.t. sphericity axes 


1.0 


In-plane p± in GeV w.r.t. thrust axes 


1.0 


Thrust minor, m 


2.0 


C parameter 


1.0 


Scaled momentum, x p — |p|/|pbeam| 


1.0 



Tab. 6: Observables used in the tuning and associated weights for observables 



taken from 45 . 



Observable 


Weight 


Observable 


Weight 


Mean p u (770) multiplicity 


10.0 


Mean Xd(3510) multiplicity 


10.0 


Mean A++(1232) multiplicity 


10.0 


Mean D + multiplicity 


10.0 


Mean K*+(&92) multiplicity 


10.0 


Mean £ + multiplicity 


10.0 


Mean £° multiplicity 


10.0 


Mean A (1285) multiplicity 


10.0 


Mean A° multiplicity 


10.0 


Mean / 2 (1270) multiplicity 


10.0 


Mean K + multiplicity 


10.0 


Mean J/ip(lS) multiplicity 


10.0 


Mean S°(1530) multiplicity 


10.0 


Mean £? + multiplicity 


10.0 


Mean A(1520) multiplicity 


10.0 


Mean B** multiplicity 


10.0 


Mean D*+(2112) multiplicity 


10.0 


Mean A+ multiplicity 


10.0 


Mean £"(1385) multiplicity 


10.0 


Mean D° multiplicity 


10.0 


Mean /i(1420) multiplicity 


10.0 


Mean /2(1525) multiplicity 


10.0 


Mean 0(1020) multiplicity 


10.0 


Mean T, ± multiplicity 


10.0 


Mean #|°(1430) multiplicity 


10.0 


Mean D^ 2 multiplicity 


10.0 


Mean Q~ multiplicity 


10.0 


Mean K*° (892) multiplicity 


10.0 


Mean £±(1385) multiplicity 


10.0 


Mean £~ multiplicity 


10.0 


Mean ip(2S) multiplicity 


10.0 


Mean ir + multiplicity 


10.0 


Mean L>*+(2010) multiplicity 


10.0 


Mean jo (980) multiplicity 


10.0 


Mean B* multiplicity 


10.0 


Mean £+(1385) multiplicity 


10.0 


Mean ir° multiplicity 


10.0 


Mean D + s multiplicity 


10.0 


Mean rj multiplicity 


10.0 


Mean p multiplicity 


10.0 


Mean a + (980) multiplicity 


10.0 


Mean B® multiplicity 


10.0 


Mean multiplicity 


10.0 


Mean K° multiplicity 


10.0 


Mean p+(770) multiplicity 


10.0 


Mean B + ,B^ multiplicity 


10.0 


Mean EE - multiplicity 


10.0 


Mean A multiplicity 


10.0 


Mean w(782) multiplicity 


10.0 


Mean ?/(958) multiplicity 


10.0 


Mean T(IS') multiplicity 


10.0 







Tab. 7: Multiplicities used in the tuning and associated weights for observables 
taken from 53 
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Observable 



b quark fragmentation function f(x^ cak ) 

Mean of b quark fragmentation function /(a;™ ak ) 



Weight 



7.0 
3.0 



Tab. 8: Observables used in the tuning and associated weights for observables 



taken from 54 . 



B Simulation of h° bb using the POWHEG Method 

The NLO differential decay rate in the POWHEG |39j approach is 

da = B($ m )d$ B 



A£ LO (0) + °(p T )^^d«I> 1 



(3) 



where 



B($ m ) = B($ m ) + V($ t 




R($ m , $0 - A($™ $ i) d$i. (4) 



Here _B($ m ) is the leading-order Born differential decay rate, V($ m ) the regu- 
larized virtual contribution, _Di($ m ,$i) the counter terms regularizing the real 
emission and -R($ m , $i) the real emission contribution. The leading-order process 
has m outgoing partons, with associated phase space $ m . The virtual and Born 
contributions depend only on this m-body phase space. The real emission phase 
space, $ m +i, is factorised into the m-body phase space and the phase space, $i, 
describing the radiation of an extra parton. 

The Sudakov form factor in the POWHEG method is 



A NLO 



exp 



Pt) 



(5) 



where ^($ ra , $i) is the transverse momentum of the emitted parton. 

In order to implement the decay of the Higgs boson in the POWHEG scheme 
in Herwig++ we need to generate the Born configuration according to Eq.|4]and 
the subsequent hardest emission according to Eq.[5j The generation of the trun- 
cated and vetoed parton showers from these configurations then proceeds as de- 
scribed in Refs. |§|5T||55||56| . 

The virtual contribution for h° — > bb was calculated in Ref. 57 . The corre- 
sponding real emission contribution, see Fig. M is 



\M R \ 2 = \M 2 \ 



CF87ra s 



M£(l 



4/x 2 



X Cl 



+ 



(8/i 4 - 6/i 2 + 1) 

1 -X q )(l -X q ) 



2(1 - 4// 



X n 



2/i 2 (l -4^ 2 ) 



1 



X Cl 



(6) 
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Fig. 9: The two real-emission processes contributing to the NLO decay rate. 



where A^2 is the leading-order matrix element, Cf = 3, rn q is the mass of the 



bottom quark, M H is the mass of the Higgs boson, /1 



m and x * = si- we 



use the Catani-Seymour subtraction scheme 58 where the counter terms are 
8tt at s , 



Di = C T 



■\M 2 \ 2 



x 



2(1 - 2/x 2 



Xj 1 2 



1 - xj - 2/i 2 
x) - V 1 - 2^ 2 



2 + 



.1 j 



2/i 



1 2/i 2 



.7 J 



(7) 



where for Dj, i is the emitting parton and j is the spectator parton. In practice, 
as the counter terms can become negative in some regions, we use 



El^($ m ,$i)| 



We have also regulated singularities in the virtual term V($ m ) with the inte- 
grated counter terms from the Catani-Seymour subtraction scheme allowing us 
to generate the Born configuration according to _B($ m ). 

The hardest emission for each leg is generated according to 



a nlo _ 

^iR — ex P 



M 2 H 



16tt 2 (1 - 4// 
dxi dx 2 d(f) 



■x 





Di 






j 


D 3 \ 



0(hr($ m ,$i)-PT) 



(9) 



However this form is not suitable for the generation of the hardest emission. 
Instead we perform a Jacobian transformation and use the transverse momentum, 
Pt, rapidity, y, and azimuthal angle, 0, of the radiated gluon to define the phase 
space $1. 
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The momenta of the Higgs boson decay products are 



Pi 

P2 
P3 



2 

Mh 
2 



x\, — x± cos(0), — x± sin(0), i-y/a^" 
x 2 ;0,0,- 



X i 



4/r 



x 3 ; x± cos(4>),x± sm(4>), ±\J x"l — x 



i / 



(10a) 
(10b) 
(10c) 



where partons 1,2,3 are the radiating bottom quark, spectator antibottom quark 
and radiated gluon, respectively. The energy fractions Xi = and 
x ± — ff^- Using the conservation of momentum in the ^-direction and 
xi + X2 + X3 = 2 gives 

2\2 



2 / n x2 (-2 + 2xi + 2x 2 - x 2 x x - x 2 " 
x ± = (2 - xi - x 2 ) 



x 2 2 — 4/i 2 

Together with the definition, x 3 = x± cosh y, we obtain the Jacobian 



dxidx^ 



x ± 



0C j_ (yX ^ 



4/x 2 



'ir. 



(12) 



dp T dy M H (xix 2 - 2fi 2 (xi + x 2 ) + x 2 - x 2 ) ' 

for the transformation of the radiation variables. 

We can then generate the additional radiation according to Eq.|9] using the 
veto algorithm [3]. To achieve this we use an overestimate of the integrand in 



the Sudakov form factor, f{pr) = — , where c is a suitable constant. We first 



generate an emission according to 
ATM = exp 
using this overestimate, where y m; 



cosh" 



J/max 



dp T dy- 



Pt 



1 Mi 



2p™ m 



, Vn 



(13) 



p™ ax is the 



maximum possible transverse momentum of the gluon and p™ m is a parameter 
set in the model, taken to be 1 GeV. 

The trial value of the transverse momentum is obtained by solving 1Z = A^ or , 
where TZ is a random number in [0, 1], i.e. 



p T = p™ X 7£c(2/m ax -y min ) 



(14) 



Once the trial px has been generated, y and are also generated uniformly 
between [y m i n , y ma x] an d [0,27r], respectively. The energy fractions of the partons 
are obtained using the definition x 3 = x± coshy, 



2(^3 - I) 



r 2 
2 



X 1 



3a;3 - 2 + -^x 3 - x\ - x" 3 



± - 4X0*3 - 1)(V + x 3 -l)- /i 2 xi)|. 



(15) 
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and X2 using energy conservation. As there are two solutions for x\ both solutions 
must be kept and used to calculate the weight for a particular trial px- The signs 
of the z-components of the momenta are fixed by the sign of the rapidity and 
momentum conservation. Any momentum configurations outside of the physically 
allowed phase space are rejected and a new set of variables generated. The 
momentum configuration is accepted with a probability given by the ratio of 
the true integrand to the overestimated value. If the configuration is rejected, 
the procedure continues with p™ ax set to the rejected px until the trial value of 
Pt is accepted or falls below the minimum allowed value, p™ m . This procedure 
generates the radiation variables correctly as shown in Ref. [3]. 

This procedure is used to generate a trial emission from both the bottom 
and antibottom. The hardest potential emission is then selected which correctly 
generates events according to Eq. [9] using this competition algorithm. 
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